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Abstract 

We study the minimum-time orbit phasing maneuver problem for a constant-current electrodynamic 
tether (EDT). The EDT is assumed to be a point mass and the electromagnetic forces acting on the 
tether are always perpendicular to the local magnetic field. After deriving and non-dimensionalizing the 
equations of motion, the only input parameters become current and the phase angle. Solution examples, 
including initial Lagrange costates, time of flight, thrust plots, and thrust angle profiles, are given for a 
wide range of current magnitudes and phase angles. The two-dimensional cases presented use a non-tilted 
magnetic dipole model, and the solutions are compared to existing literature. We are able to compare 
similar trajectories for a constant thrust phasing maneuver and we find that the time of flight is longer 
for the constant thrust case with similar initial thrust values and phase angles. Full three-dimensional 
solutions, which use a titled magnetic dipole model, are also analyzed for orbits with small inclinations. 
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Reference orbit radius 
Magnetic Field Vector 
Derivative of State Vector 
Force vector 
Magnetic Field Frame 
Hamiltonian 
Total dipole strength 
Current 

Tether direction vector 
Mass 

Dipole direction unit vector 
Reference mean orbit velocity 
Inertial Coordinate Frame 
Orbit radius 


R 0 Radius of the earth 
R Position unit vector 
T Dimensionless thrust 

Tq Dimensionless thrust at original orbit radius 
TU Time Unit 
v x ,v y , v z Velocity components 
x,y,z Position components 
/3 Thrust Coefficient 
Phase angle 
Aq Lagrange costate 

fi Gravitational parameter for attractive body 
ijj Thrust control angle 
p Lagrange costate ratio 
() Inertial time derivative 


1 Introduction 

Using an electrodynamic tether (EDT) as a means of propulsion is a relatively new concept that has not 
been fully developed. There is much in the literature concerning the trajectories of EDTs, 15 and Williams 
gives some solutions for optimal control of an EDT. 1 However, the equations of motion developed in the 
literature are based on dimensional parameters, making the results specific to the cases presented. Also, 
in many studies involving EDTs, 1,4 the magnetic field model is a non-tilted dipole. There are no other 
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known studies that use a tilted dipole model in an optimal control problem. Some other assumptions include 
restricting the EDT orbit to the equatorial plane and aligning the EDT with the nadir, or radial, direction. 
This is an unrealistic assumption for any practical application. 2-5 All of these assumptions do not produce 
enough information create a usable control model. 


2 Model 

This section will describe the basic physics behind the EDT and how the physics leads to the development 
of the coordinate system used in analysis. 

2.1 Physics 

The EDT creates thrust by using a Lorentz force, which is defined by the following equation, 


F = iLxB (1) 

where F is the force, i a scalar current value, L is the tether vector, and B is the magnetic field vector. As the 
conductive tether moves through the earth’s magnetic field, the field induces a current through the wire to 
counter the change in magnetic flux passing through the tether. The induced current creates a Lorentz force 
in the opposite direction of flight, slowing the EDT and reducing its orbit altitude. 6 Essentially the tether is 
converting orbit energy to electrical power in the tether. If the EDT supplies its own current, it can create 
a force in the direction of travel to accelerate the EDT and raise its altitude. Because the tether is not a 
circuit, it cannot complete the circuit on its own. It needs free electrons from the earth’s ionosphere. Devices 
on each end of the tether exchange electrons from the ambient plasma. An illustration of this process is seen 
in Figure 1. Because the EDT must operate in the ionosphere, the orbit altitude must be below about 2000 
km. However, this range includes all current manned missions and all other space applications requiring a 
low earth orbit. 



Figure 1: Physics Behind EDT Operation 7 


2.2 Coordinate System and Magnetic Field Model 

This optimization problem uses two coordinate frames: an Earth-centered inertial frame, or ECI frame, and 
a magnetic field frame, F. The ECI frame is defined as N : hi, h 2 , n 3 , where hi is in the Vernal Equinox 
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direction, 113 is parallel to the earth’s axis of rotation, and n 2 is the cross product of 113 and ni. The 
magnetic field frame, F : 6 i, 6 2 ,o 3 is described below. 

The B vector points in the direction of the earth’s magnetic field at a given point in space, as shown in 
Figure 2. The magnetic field used in this problem is based on the tilted dipole model. 8 In vector form, the 
direction and magnitude of the magnetic field are given by 

B(R) = (3 (m • R) R - m) (2) 

where R 0 is the radius of the earth, Ho is the total dipole strength, m is the dipole direction, and R is 
the normalized position vector. The International Association of Geomagnetism and Aeronomy defines the 
dipole direction and total magnetic field strength in the International Geomagnetic Reference Field report. 8 
The unit vector of B defines 63 . 
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Figure 2: Tilted Dipole Magnetic Field Model 

The two remaining unit vectors of the magnetic field frame, 61 and o 2 , require two reference planes to 
define them. The first plane contains both rii and the dipole direction vector, m. The second plane is 
perpendicular to 63. The line of nodes created by the intersection of these two planes defines 61, as depicted 
in Figure 3. The mathematical definition of 61 is 


61 = 63 x (m x fi i) ( 3 ) 

The vector o 2 is defined by o 2 = 63 x 61. 

The thrust control angle, ip, gives the direction of the force in the 61, o 2 plane. Solving for the control 
angle history as a function of time is the main goal of this problem. 
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Figure 3: Magnetic Field Frame Definition 
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Figure 4: Control Angle 
Definition 


2.3 Equations of Motion 

The derivation of the twelve equations of motion follows the methods given by Thorne and Hall. 9,10 We 
assume the EDT and the attracting body are point masses and that the only other external force is the 
Lorentz force. First, the second-order equations of motion are written as a system of first-order differential 
equations 
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where r = x 2 + y 2 + z 2 , /q is the scalar component of a perturbation force in the r, y, or z direction, 
and /i is the gravitational parameter for the attracting body. Here, x, y , and z are components of the EDT 
position vector with respect to the origin of the ECI frame, TV, and v x , v y , and v z are the derivatives of x, y , 
and z with respect to the inertial frame. To simplify the analysis, the magnetic force vector is represented 
in F frame coordinates. A coordinate transformation is required to rotate this vector to the inertial frame. 
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where Cij is an element in the rotation matrix, [R] . These force components are inserted into Equations 7-9 
to form 


v x = - H — (Cu cos ip + C 12 sin VO 

r 6 m 
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Equations 4-6 and 13-15 are then used to define the Hamiltonian, H = A T d, where A is the Lagrange 
costate vector defined as [\ x , A y , A^, A Vx , X Vy , X Vz ] T and d is a state vector containing the components of the 
derivatives of position and velocity. The control law is then obtained by taking the partial derivative of H 
with respect to ip and setting it equal to zero 


dH 

= \x (— Cu sin V^ + Ci 2 cos VO + K y (— C 21 sin -0 + C 22 cos ip) + X vz (— C 3i sin^ + C 32 cos VO = 0 (16) 


which leads to the optimal thrust angle in terms of the costates 
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This control law is inserted into Equations 13-15 to get the components of velocity in terms of the Lagrange 
costates: 
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To complete the analysis, we need six more equations to solve a system with twelve state variables. The 
remaining equations of motion are 
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( 25 ) 
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We wish to use non-dimensional units in the analysis so that the results are universal and so that the time 
of flight is simply the change in true anomaly during the maneuver. To complete the transformation, /i is 
set to unity, and BiL/m is represented by a dimensionless thrust variable, T, defined as 


BiL/m 
mR 0 n 2 


( 28 ) 


where n is the dimensional circular orbit velocity at R Q . 

The equations of motion are constructed so that they describe minimum-time motion. In other words, 
the EDT, whose trajectory is described by the equations of motion, arrives at any given position in the 
minimum time possible. This formulation makes time of flight a cost function in the optimization problem. 
The other cost function used is the difference between the position and velocity components of the EDT at 
t = tf and the target position and velocity components at t = tf. The other cost function is proportional to 
how well the EDT matches the target spacecraft’s position and velocity. The difference between the position 
and velocity of the EDT and that of the target spacecraft is called the residual. The EDT meets its target 
when the residual approaches zero. 


2.4 Algorithm 

The first step of the algorithm uses a numerical method called simulated annealing (SA). It mimics the 
physical process of taking a highly unorganized, melted material and reducing the temperature slowly so 
that the crystalline structure is nearly perfect. Simulated annealing involves lowering the energy state, or 
performance index, of the system slowly so that the algorithm can find possible minima. As the energy state 
lowers, the number of possible locations for the global minima within the search space decreases until only 
one location remains. The exact algorithm for simulated annealing can be found in Refs . 11 and . 12 

Simulated annealing is used to find initial costates and times of flight that produce the smallest residual. 
The problem with SA is that it has a low convergence rate and low accuracy. The accuracy is refined 
by using a root-finder, such as f solve in MATLAB. Root-finders have small convergence radii but high 
convergence rates and high accuracy. Finally, once one solution is found, a family of solutions can be found 
using Newton’s method of continuation . 11 

To obtain full three-dimensional optimal control solutions, comparable two-dimensional solutions must 
be found first. In the two-dimensional problem, the target orbit is in the earth equatorial plane and m points 
in the — h 3 direction. Having the dipole in the — h 3 direction keeps the electromagnetic force in the hi — h 2 
plane because B is always perpendicular to the plane. Also in the two-dimensional case, the N frame and 
the F frame have the same orientation because of how each frame is defined. 

The SA algorithm finds an initial guess for a two-dimensional system, which is then refined to a solution 
within a specified tolerance. The two-dimensional solution is used as an initial guess for a full tilted dipole 
problem in three-dimensional space. 


3 Results 

To gain insight into the EDT optimal control problem, we start with the simple two-dimensional case and 
compare it to the literature. In the two-dimensional case, the magnetic field model used is a non-tilted 
dipole, and all orbits are in the equatorial plane. We then use the two-dimensional solutions to develop 
three-dimensional trajectories. 
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3.1 Two-Dimensional Optimal Control Problem 

Several conclusions can be drawn from studying the two-dimensional problem. As a reference, the following 
results are compared to the work of Hall and Collazo-Perez. 9 Their work focused on constant thrust, 
minimum-time optimal phasing maneuvers in two-dimensional space. The authors provided several solutions 
that were representative of their work and trends that encompassed their results. Two aspects of their work 
can be directly compared to the two-dimensional EDT optimal control problem. One comparison is obviously 
the trajectory. The other is the near- invariance of the initial costates when changing the dimensionless thrust 
and phase angle by an order of magnitude. 

3.1.1 Trajectory Analysis 

When comparing the optimal trajectory of a constant current EDT with that of a constant thrust spacecraft, 
intuition suggests that the EDT will catch up to its target faster because thrust increases when the orbit 
radius decreases, in accordance with (1). A higher thrust should propel the EDT to the target more quickly 
than a constant thrust model. 

Calculations show that the EDT does indeed arrive at its target more quickly. Figures 5 and 6 illustrate 
a comparison between the trajectories for an EDT and a spacecraft with constant thrust. The red, dotted 
line is the original orbit and the blue, solid line is the trajectory. In the case depicted, T 0 = 0.5 and </> = 1.46 
rad, where To is the dimensionless thrust at the original orbit radius. The value for thrust is unusually large 
for an EDT, but this example is chosen so that differences in their trajectories would be obvious on a plot 
of the transfer trajectories. The plot shows the EDT arriving at its target about 0.3 TU before the constant 
thrust spacecraft arrives. Also the shape of the path at mid- flight is more rounded for the constant thrust 
trajectory than for the EDT trajectory. As stated by Equations 1-2, thrust is proportional to R~ 3 . The 
increase in thrust give the EDT the ability to fly almost directly at its target, thus creating the flattened 
shape in its trajectory. 




Figure 5: EDT Trajectory Figure 6: Constant Thrust Trajectory 

As stated before, a large phase angle and a large thrust illustrate the EDT’s shorter time of flight, but 
the time of flight for the EDT is shorter than the constant thrust spacecraft for all cases. When comparing 
initial costates presented in Ref. 9 with the EDT initial costates for several cases, the EDT reaches its target 
faster in every case. Table 2 illustrates this fact. The difference in time of flight is more pronounced in large 
phase angle maneuvers. For the case with = 0.54, the difference in time of flight is 0.0882 Time Units (TU) 
while the difference in time of flight is only 0.001525 TU for = 0.0074 rad. The small phase angle does 
not allow the EDT to dip far enough into the orbit to use its advantage over the constant thrust spacecraft. 

Figures 7 and 8 show the control angle, / 0, for each spacecraft. For the two-dimensional case, -0 is 
measured with respect to hi. The shapes are similar with the control angles starting in the third quadrant 
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Table 2: Lagrange Costate Comparison 




IO 

o 

II 

E? 

IO 

o 

II 

E? 

T 0 = 0.05, 

T 0 = 0.005, 

T 0 = 0.005, 

T = 0.005 



(j) = 1.46 

4> = 0.54 

(j> = 0.89 

4> = 0.0074 

0 = 0.022 

(j) = 0.1 

EDT 

^ yo 

0.42896 

0.099345 

0.34992 

0.18569 

0.43658 

0.33357 


^v x0 

0.62643 

0.65774 

0.44331 

0.72805 

0.69926 

0.43892 


^ VyO 

0.42613 

-0.039660 

0.98215 

0.000697 

0.49590 

0.99619 


tf 

2.51007 

1.95232 

5.42443 

2.45092 

3.78494 

6.16831 

Constant 

^ yo 

0.47313 

0.14480 

0.34156 

0.18627 

0.43713 

0.33270 

Thrust 

^ VxO 

0.63251 

0.66574 

0.43092 

0.72812 

0.69905 

0.43752 


^ VyO 

0.49960 

0.00012 

1.00027 

0.001362 

0.49755 

0.99824 



2.78398 

2.04048 

5.58198 

2.45245 

3.79228 

6.18639 


to lower the orbit and increase speed. The thrust then rotates 180° at the half way point to raise the orbit 
to meet with the target. The control angle rotates in this manner for all cases studied. 




Figure 7: EDT Control Angle Figure 8: Constant Thrust Control Angle 


3.1.2 Near-Invariance of the Initial Lagrange Costates 

One important discovery in Ref. 9 is that when T and (j) are both increased by an order of magnitude, the 
values of the initial Lagrange costates and the time of flight change by a small amount for moderate and low 
thrust cases. This near-invariance of the initial costates is important because it allows for the calculation of 
a whole family of solutions once one solution is found. This trait is also present in the EDT two-dimensional 
optimal control problem. 

Table 3 gives three representative examples that validate the claim of near- invariance for moderate and 
low thrust cases and shows its limitations. The first example, the case with the lowest thrust value in the 
group, has the least amount of variance as T and </> are increased by an order of magnitude. The costates 
and the time of flight are constant to three decimal places. The second example shows an invariance between 
the two costate sets with a precision of two significant digits. The last example compares a high thrust, T 
= 0.5, with a thrust one order of magnitude smaller. This comparison shows the most variation of all the 
three cases. Further exploration shows large variations in the initial costates when large values of T and 
are decreased by an order of magnitude. The large variations in the costates are caused by the different 
ranges of thrust magnitudes experienced during the maneuver. There is a larger difference in the thrust 
range experienced during the maneuver in Example 3 than there is during the maneuver in Example 1. A 
larger thrust range causes a greater variance between costate sets. 
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Table 3: Demonstration of Near-Invariance of the Initial Costates 



T 

0 

^ yo 

^ V X 0 

^ VyO 


1 

0.0005 

0.00074 

0.18527 

0.72864 

-0.0003115 

2.45355 


0.005 

0.0074 

0.18569 

0.72805 

0.000697 

2.45092 

2 

0.009197 

0.001 

-2.52311 

0.34754 

-0.85886 

0.66334 


0.09197 

0.01 

-2.48727 

0.34688 

-0.84637 

0.66314 

3 

0.05 

0.146 

0.38861 

0.71940 

0.32728 

3.18374 


0.5 

1.46 

0.42896 

0.62643 

0.42613 

2.51007 


3.2 Three-Dimensional Optimal Control Problem 

Now that two-dimensional solutions for the EDT optimal control problem have been established, these 
solutions can be used as initial guesses to full three-dimensional solutions using the tilted-dipole model. 
These initial guesses prove to only be useful to find zero inclination phasing maneuvers. Also not all two- 
dimensional solutions can be used as initial guesses. Solutions with phase angles over 0.01 radians or thrusts 
over 0.01 do not converge when using f solve. When a three-dimensional, zero inclination angle solution is 
found, the reason behind the small convergence radius becomes clear. The F frame rotates greatly with time 
relative to the inertial frame due to the changes in B, and by extension the o 3 vector. Because the thrust 
vector is constrained to 6i — o 2 plane, which can rotate rapidly with the passage of time, the algorithm 
has difficulty converging to an solution. Once a solution is found, one must consider the dynamics of the 
solution. Because the thrust is always perpendicular to the magnetic field, the EDT generates the most 
thrust when it is in the 6i — o 2 plane. For this optimal control to work, the tether must not only rotate 
about o 3 , but it must also rotate as the F frame rotates. 

3.3 Case Study: Changing Orbit Plane Inclination Angle 

In this section, we analyze the changes of the initial costates and the time of flight as the orbit inclination 
increases from 0°-7.5°. The constant variables in this case study are the phase angle and the product of 
the tether current and tether length, iL. The values of these variables are 0.002 radians and 300 A-m, 
respectively. The product of the current and tether length are non-dimensionalized using (28). 

The task of finding phasing maneuver solutions for increasing inclination angles is a difficult one because 
of the trajectory’s high sensitivity to the initial guess. Changing an initial guess by 1% can cause the tether 
to be off target by a distance on the order of 100 km. The high sensitivity is due to the changing magnetic 
environment caused by the tilted magnetic dipole. 

The sensitive nature to the initial guess is not reflected in all initial conditions. The appearance of the 
plots for each initial condition as a function of the inclination angle can be divided into two groups. Figure 
3.3 represents Group 1, which includes X y and X Vx , and Figure 3.3 represents Group 2, which includes A^, 
X Vy , X Vz costates and tf. The slopes of the data in Group 1 change rapidly as i increases, making the 
task of predicting new solutions using continuation difficult. On the other hand, the new solutions for the 
initial conditions in Group 2 are easy to predict because the slope changes gradually. The data in Group 2 
can be curve-fitted using a linear regression technique and extruded to make a new prediction for a higher 
inclination angle. 

3.3.1 Control Angle Analysis 

Two distinct control angle patterns appear for the inclination angle range analyzed in this paper. For 
inclination angles between 0 and 2 degrees, the control angle plot looks like Figure 11. The overall trend is 
that the angle continuously decreases over time. A similar trend is seen for the half-orbit 2D solutions as 
presented in Section 3.1.1, which makes sense because having inclination angles below 2° creates a similar 
magnetic environment to that in the 2D problem. 

For inclination angles between 3 and 7.5 degrees, the control angle plot looks like Figure 12. Here ^ is 
in the third quadrant for almost half the manuever in order get inside the target orbit to catch up. The 
control angle then spins around with an increasing ^ to the first quadrant to slow down before finishing the 
maneuver in the fourth quadrant. If the two control angle plots are compared, the figures show that the 
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Figure 9: Variations of Group 1 with Increasing i Figure 10: Variations of Group 2 with Increasing i 



Figure 11: Control Angle for i = 1° Figure 12: Control Angle for i = 5° 


control angles beyond 1.8 TU are almost identical. The only difference in 'ip is the direction that the EDT 
spins after the initial portion of the maneuver. This change in spin direction is also apparent in Ref. 9 when 
a parameter, in that case an initial costate, is changed by a small amount. 


4 Conclusion 

Solutions to the minimum-time, constant-current orbit phasing maneuver for an electrodynamic tether have 
been found for a wide range of two-dimensional cases as well as limited three-dimensional cases. The 
two-dimensional results are compared with existing data from the literature to provide insight into the 
characteristics of the solutions. In all cases, the EDT reaches its target faster than a constant thrust 
spacecraft because the Lorentz force acting on the tether increases as orbit radius decreases. Also sets of 
initial Lagrange costates, whose parameters T and (p each vary by an order of magnitude, are shown to be 
near- invariant for moderate to low values of thrust and phase angle. These two-dimensional solutions are 
used as initial conditions to solve for three-dimensional solutions which include a tilted-dipole magnetic field 
model. These initial conditions are only useful for zero inclination orbits. The most inclined orbit phasing 
maneuver solved at the time of publication is at 5°. The reason for the difficulty in converging to a more 
inclined solution is due to the rotation of the F frame. 
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Further research is needed to fully develop three-dimensional solutions for inclination angles over 7.5° 
and for phase angles and current values different than those analyzed in this paper. The behavior of the 
F frame becomes more erratic as inclination angle increases and must be taken into account when finding 
additional solutions. After more solutions have been found, the same solution algorithm presented can be 
used to find optimal trajectories for other combinations of and iL and for other simple maneuvers such 
as inclination angle change and orbit radius change. These basic maneuvers must be solved in order for 
electrodynamic tethers to become practical. 
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